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Abstract. The following work compares two popular mixed finite elements used to model subsur- 
face flow and transport in heterogeneous porous media; the lowest order Raviart-Thomas element 
and the variational multiscale stabilized element. Comparison is made based on performance for 
several problems of engineering relevance that involve highly heterogenous material properties (per- 
meability ratios of up to 1x10^), open flow boundary conditions (pressure driven flows), and large 
scale domains in two dimensions. Numerical experiments are performed to show the degree to which 
mass conservation is violated when a flow field computed using either element is used as the advec- 
tion velocity in a transport model. The results reveal that the variational multiscale element shows 
considerable mass production or loss for problems that involve flow tangential to layers of differing 
permeability, but marginal violation of local mass balance for problems of less orthogonality in the 
permeability. The results are useful in establishing rudimentary estimates of the error produced by 
using the variational mutliscale element for several different types of problems. 



1. INTRODUCTION 

In the field of subsurface modefing, choosing the right numerical formulation for a particular 
problem of interest can be a difficult choice. There are a number of trade-offs that one must consider. 
Two very important considerations involve accuracy and computational efficiency. For problems 
that involve species transport (reactive or not) through a porous domain, artificial production or loss 
of constituents can be detrimental, especially in the case of transport of contaminants or pollutants. 
At the same time, accuracy comes with a cost, which manifests in the form of computation time, 
mesh resolution requirements, and/or time step restrictions. An astute analyst must balance these 
tradeoffs such that sufficiently accurate results are obtained at a reasonable cost. 

In this paper, we compare two popular mixed finite element formulations to solve the fiow equa- 
tions in heterogeneous porous media. The two mixed formulations are the lowest order Raviart- 
Thomas formulation (RTO) [IJ, and the variational multiscale stabilized formulation (VMS), which 
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has been studied in references [SUHKHIS]. Mixed methods, that treat velocity and pressure as inde- 
pendent variables, as opposed to single field formulations, provide a more accurate representation 
of the velocity field, but are susceptible to instability engendered by the saddle-point nature of the 
system of equations [6j. The VMS element evaluated in this work consists of equal order linear 
interpolation of the velocity and pressure fields, for which the degrees of freedom live at the node 
locations. The VMS element achieves stability through the addition of residual based terms to 
the formulation such that the LBB condition is circumvented. The lowest order Raviart-Thomas 
finite element formulation (RTO) consists of evaluating fluxes over the midpoint of element edges 
and constant pressures within the element and is based on an LBB stable finite element space. 
A schematic of both elements is presented in Figure [T} Both elements are used to solve several 
realistic engineering problems involving flow and transport in heterogeneous porous media. 

The pursuit of stable, locally conservative finite element formulations has received substantial 
attention (for examples see [TJ [U [9] among many others). In addition, a number of comparisons 
have been made between various formulations (for example, see \10 \ IT U IT2]). Of particular interest 
is the work presented in |12j that showed a computational efficiency advantage for control volume 
approaches vs. mixed methods for problems with only moderate heterogeneity. 

It is well known that the VMS element, which is based on a Galerkin formulation is not locally 
mass conserving (without additional sophistication being built into the algorithm as in |13| [TT] ) , 
but the degree to which a realistic simulation will suffer due to this shortcoming of the VMS 
element is not well documented. Additionally, although the RTO element is locally conservative, 
it has only been shown viable for triangular or tetrahedral elements (quadrilateral variants of this 
method do not pass the patch test). In the context of heterogeneous materials, this work aims, to 
evaluate the strengths and weaknesses of either approach for problems with highly varying material 
permeability and predominantly open flow boundary conditions. The results have applications in 
subsurface simulations for which the far-field boundary is typically described with some type of 
pressure boundary condition. 

1.1. Main contributions of this paper. One of the main goals of this paper is to assess the per- 
formance of the variational multiscale formulation with respect to local (or element) mass balance 
using a coupled flow-transport model. The assessment is performed by comparing the variational 
multiscale formulation with the lowest-order Raviart-Thomas formulation, which possesses the lo- 
cal mass balance property. We also take into account the dependence of viscosity on the pressure, 
which is the case with most organic liquids and supercritical carbon-dioxide. Using Barus model 
(in which the viscosity varies exponentially with respect to pressure), we study the affect of the 
dependence of viscosity of pressure on the error in the local mass balance under the variational 
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multiscale formulation. Using canonical problems, we also compare the performance of the vari- 
ational multiscale formulation and the lowest order Raviart-Thomas formulation with respect to 
memory storage, time for assembling matrices, and total CPU time. 



1.2. Organization of the paper. The remainder of this paper is organized as follows. In Sec- 
tion [2j we present a coupled flow-transport mathematical model. In Section |3j we present weak 
formulations to solve both flow and transport problems. For the flow problem, we consider the 
variational multiscale formulation and the lowest-order Raviart-Thomas formulation. In Section |4| 
we present numerical results for flow and transport for three canonical problems: a stratified set 
of layers with differing permeabilities, a cylindrical inclusion of high permeability in a domain of 
lower permeability, and a leaky well scenario where a storage aquifer leaks the injected constituent 
through an abandoned well. In Section [5| we conclude with a summary of the important research 
findings of this paper. 

2. A COUPLED FLOW-TRANSPORT MATHEMATICAL MODEL 

Consider the transport of a chemical species in a rigid porous medium (that is, the deformation 
of the porous solid is neglected). We shall assume that the porous medium is fully saturated with 
a (moving) fluid. Herein, we consider modified Darcy equation for the flow problem which (as the 
name suggests) is a modification to the standard Darcy equation [TS] . This model takes into account 
the dependence of viscosity on the pressure. The fluid is still assumed to be incompressible. The 
transport model consists of transient advection-diffusion equation wherein the advection velocity 
is given by the flow problem. 

2.1. Governing equations for the flow problem. Let O be a bounded open domain, and F 
be its piecewise smooth boundary (F := — i}, where 0, is the set closure of f]). A spatial point 
is denoted hy x £ ^1. The gradient and divergence operators with respect to x are, respectively, 
denoted by grad[-] and div[-]. Let p : — )• M denote the pressure, and the velocity field is denoted 
by i» : — )• M"*^, where "nd" is the number of spatial dimensions. We denote the density by p, the 
permeability by k{x), and the coefficient of viscosity by p. We consider the following expressions 
for the dependence of viscosity on the pressure: 

fi{p) = po (la) 
fi{p) = fioexp[j3p] (lb) 

where //q > is the coefficient of viscosity at a reference pressure (e.g., at the atmospheric pressure). 



and /3 > is a constant. Equation (la) assumes the viscosity is constant, which is the case in the 



standard Darcy fiow model fT5] . Equation ( lb ) is commonly referred to as Barus formula (161 
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We write the governing equations for the flow in first-order form (or mixed form) in which the 
velocity and pressure are treated as independent variables. The governing equations for the flow 
problem can be written as follows: 

in n (2a) 

in n (2b) 

on r (2c) 

where b{x) is the specific body force, (j){x) is the prescribed volumetric flow rate, il){x) is the 
prescribed normal component of the velocity on the boundary, and n{x) is the unit outward normal 
vector to T. Note that for a given ft (and hence for a given F), (^(cc) and ip{x) are not independent 
of each other, and they need to satisfy the following compatibility condition: 

/ (l){x) dn= [ ip{x) dr (3) 
Jn Jt 

Since we have prescribed only the normal component of the velocity on the whole boundary, one 
need to enforce one more condition to obtain unique solution in the pressure field. This is typically 
achieved by meeting the following condition: 

/ p{x) dn = (4) 
Jn 

Another way to achieve uniqueness in the pressure field is to prescribe the pressure at a point 
to fix the datum, which is computationally the most convenient and is employed in this paper. 
In pressure-driven flows, there is no need to enforce additional conditions for uniqueness as the 
pressure is prescribed on a set of non-zero measure on the boundary. 

2.2. Governing equations for the transport problem. Let c{x, t) denote the concentration, 
and X denote the time interval of interest. The boundary T is divided into two complementary 
parts and F^. F^ is that part of the boundary on which Dirichlet boundary condition is 
prescribed (that is, the concentration is prescribed), and F^ is that part of the boundary on which 
the Neumann boundary condition is prescribed. The transport problem is to find c{x, t) by solving 



the following initial boundary value problem: 
dc 

— + v{x) ■ grad[c] — div [D grad[c]] = f{x, t) in Q x I (5a) 

c{x,t) = cP{x,t) on F^xX (5b) 

D grad[c] • n{x) = tP{x,t) on F^xX (5c) 



where D is the constant diffusivity of the species, d^{x,t) is the prescribed Dirichlet boundary 
condition on concentration, tP(a;,t) is the prescribed flux, f(x,t) is the volumetric source term, 
and n{x) is the unit outward normal vector to F. 
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V + grad[p] = pb{x) 
divfv] = 4){x) 
v{x) ■ n{x) = tp{x) 



3. WEAK FORMULATIONS 

In this section, we present two mixed formulations for the flow problem, and a single-field formula- 
tion for solving the transport problem. The first mixed formulation is based on the Raviart-Thomas 
spaces, and the other is the variational multiscale formulation presented in Reference [5]. 

3.1. Weak mixed formulations for the flow problem. Let the weighting function correspond- 
ing to the velocity be denoted by w{x), and the weighting function corresponding to the pressure 
be denoted by q{x). The relevant function spaces for the velocity and the corresponding weighting 
function are, respectively, defined as follows: 

V := {v{x)\v{x) G {L2{n)T'^ , divH G L2 (17), trace [t; • n] = iP{x) G H-^^^{T)} (6a) 
W := {w{x)\w{x) G (L2(Jl))"'^,div[w] G L2 (17), trace • n] = on F} (6b) 

where -^2(17) is the space of square integrable functions, trace[-] is the standard trace operator 
used in functional analysis [17], and H^^^'^(T) is the dual space of H^^'^(T). In the classical mixed 
formulation, the function space for the pressure and its corresponding weighting function is defined 
as follows: 

V := {p{xMx) G L2{n), [ p{x) dn = 0} (7) 

Jn 

In the variational multiscale formulation, the function space for the pressure and its corresponding 
weighting function is defined as follows: 

V := {p{x)\p{x) G H\n), [ p{x) dO = 0} (8) 

Jo, 

where H^{Q) is a standard Sobolev space [18]. We shall employ the following notation for the 
standard inner product: 

(a;6) = / a-6dO (9) 

Jo. 

3.1.1. Classical mixed formulation and Raviart-Thomas spaces. The classical mixed formulation 
reads: Find v{x) G V and p{x) G V such that we have 

{w] j_v) - (divH;p) - {q; divH - 4>) = {w; pb) \/w{x) G W, q{x) G V (10) 

Let V/i, Wh and Vh are, respectively, appropriate subspaces of V, W and V. A (conforming) finite 
element formulation reads: Find Vh{x) G Vh and Phix) G Vh such that we have 

{wh; j^Vh) - {div[vh];Ph) - {qh;div[vh] - 0) = {wh^ph) \/wh{x) G W/i, qh{x) G Vh (H) 

It is well-known that inclusions Vh C V, W/i C W and Vh <ZV do not insure stable numerical results. 

Specifically, not all combinations of interpolation functions for the velocity and pressure are stable. 
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A mathematical theory concerning the stabiUty of formulations and choice of appropriate finite 
dimensional spaces is given by Ladyzhenskaya-Babuska-Brezzi (LBB) stability condition [iniEnillH]. 

In the literature one can find two different ways of obtaining stable results: circumventing the 
LBB condition and satisfying the LBB condition [21j. In the former category, stabilization terms 
are added to the classical mixed formulation to achieve stability. Many stabilized formulations fall 
in the former category (e.g., references [22l |23l [5]). Some of the approaches in the later category 
are Raviart-Thomas spaces [1], Brezzi-Douglas-Marini spaces |24j . In this paper, for the flow 
problem, we shall consider the lowest order Raviart-Thomas spaces and the variational multiscale 
formulation. 

We now present the lowest-order Raviart-Thomas triangular spaces. Let Th be triangulation on 
17. The lowest order finite dimensional subspaces on triangles are defined as follows: 

Vh ■■= {v = I t-J^) = aK + bxx, = ck + bKV, aK, bx, ck£R;K£ %} (12a) 

Wh ■■= {w = {w^^\w^'^^) I w'-j^'^ =aK + bxx, w^^^ = ck + bxy, ax, V, CKeR;Ke %} (12b) 

Vh ■= {pix) I p{x) = a constant on each triangle K G Th} (12c) 

3.1.2. Variational multiscale formulation. A stabilized formulation based on the variational multi- 
scale formalism can be written as follows [4J: Find v{x) G V and p{x) G Q such that we have 

{w; j^v) - {div[w];p) - {q;div[v] -(f))- {w; pb) 

- ^ ( T'^ + g^^ad[g];-(^t; + grad[p] -p6) ) =0 yw{x) £ W, q{x) £ Q (13) 
2 \k p K J 

Remark 1. If the coefficient of viscosity is constant the formulation presented in Reference [3] is 

same as the formulation discussed in References \12\ [23] . 



3.1.3. Pressure boundary conditions. It is well-known that strongly enforced outflow boundary con- 
ditions can produce spurious oscillations in the solution [251 ES EH EE]- Let Tp denote that part 
of the boundary on which the pressure is prescribed, and let r„ be that part of the boundary on 
which the normal component of the velocity is prescribed. (As usual, for well-posedness we have 
r„ U Fp = r and F„ n Fp = 0.) A simple procedure for enforcing the pressure boundary condition 
is to replace the term — (div[i(;];p) with the terms 

-(div[w;p) + {w ■ n-pQ)rp 



in the mathematical statements (10) and ([13]) and subsequent steps. In such cases, the function 
spaces for the velocity and corresponding weighting function will be modified as follows: 

V := {v{x)\v{x) G (L2(J7))'"^,divM G Ls (17), trace [-?;• n] = V'(k) G H-^/'^{T^)} 
W := {w{x)\w{x) G {L2{VL))"-'^ ,diY[w] G ^2(17), trace [tu • n] = on F^ 



Table 1. Summary of formulation properties 





Global 


Local (element) 


Convergence 


Convergence 


Formulation 


mass balance 


mass balance 


rate velocity 


rate pressure 


RTO 


Yes 


Yes 


1.0 


1.0 


VMS 


Yes 


No 


2.0 


1.0 



3.2. Weak formulation for the transport problem. We employ the standard single-field for- 
mulation to solve the transport equations, which is based on the Galerkin principle. Let us denote 
the weighting function corresponding to the concentration as w'^{x). We shall define the following 
function spaces for the concentration and the corresponding weighting function: 

Ct := {c{x, t) I c{x, t) £ H^in), cix, t) = c^ix, t) on lJ?} (14a) 
W := {w^{x) I w^ix) £ H\n),w''ix) = on fJ?} (14b) 

The weak form can be expressed as: Find c{x, t) G Ct such that 
Be 

{w; — ) + {w"; V ■ grad[c]) + (grad[u;^]; grad[c]) = {w^; /) + {w^; t^)rn yw'^ix) G (15) 
ot 

Remark 2. In addition to the local mass balance error associated with the flow problem, there 
will also be an error associated with the transport problem. Since the focus of this work is on the 
characteristics of the flow problem we are assuming that the transport error is the same for the 
velocity field generated by both the RTO and VMS formulations. 

4. REPRESENTATIVE NUMERICAL RESULTS 

In the following section we present results from three numerical examples involving flow and 
transport in porous domains with heterogeneous material properties. For all examples, the mesh 
was constructed from linear triangular elements. Identical meshes were used for both the RTO and 
VMS elements. A direct solver (Amesos umfpack [29]) was used to solve the resulting system of 
equations in serial on a desktop machine, using an Intel Xeon CPU (2.93 GHz). For all examples, 
a constant species diffusivity, D = 0.01 was used. Also, the Backward Euler method was used for 
time integration with a time step. At = 0.01 s, and the residual tolerance was set to 1 x 10^^. 

4.1. Multilayer problem. For the multilayer problem (which is a modification of an example 

given in |30[ I31|). the unit square domain Q = (0, 1) x (0, 1), shown in Figure [2| is divided into five 

layers of thickness 0.2 units, each with a unique ratio, k/fj,. Within each layer, the ratio is constant, 

but changes as a step function at the layer boundaries. The flow is driven by a pressure difference 

prescribed on the left and right boundary of the domain. The properties used for this problem are 

defined in Table [2] The computational mesh used for the multilayer problem is shown in Figure [2] 
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The exact solution for the pressure is a Unear drop across the domain. The exact solution for the 
velocity in each layer is equal to the ratio, k/ ^. Figure [s] shows the x-velocity along the centerline of 
the domain taken at a; = 0.5. Note from the figure that the RTO element captures the exact solution 
for the velocity in each layer and the jump in velocity between layers. The VMS element is not able 
to accurately capture the jump in velocity between layers. This is due to the well-known property 
that a continuous element will enforce continuity of the velocity in both the normal and tangential 
directions, although discontinuity in the tangential direction is physically realizable. Ironically (for 
a stabilized method), in addition to smoothing the velocity in the vicinity of permeability jumps, 
the VMS element shows oscillations in the velocity within each layer. It can be shown that with 
increasing disparity between layer properties, these oscillations increase in magnitude. 

Figures |4| [Sj and [6] show the results of using the computed velocity field as the advection velocity 
for the transport problem defined by equation (5a). In this case, the problem represents transport 
of a species though a stratified system of permeabilities. Along the left edge of the domain, the 
concentration is held constant in time at 1.0. As time progresses, the species advects and diffuses 
though the domain at different rates, given the local flow velocity with each layer. Even though 
the VMS element shows significant overshoots and undershoots for the velocity in the vicinity of 
the permeability jumps, the two elements show close correlation in terms of concentration of the 
species in time. Figure |4] shows that the VMS element exhibits latency in tracking the species front 
in time, but Figure [6] suggests that the overall concentration, or the integral of c over the domain 
is preserved. There is only a marginal difference between the two elements in predicting the time 
at which a steady state solution is reached. Also, both elements capture a smooth variation of 
concentration along the front, even thought the VMS velocity field is not smooth. 

Figure [7] shows the degree to which the VMS element violates local mass balance for problems of 
this nature. (A similar figure is not presented for the RTO element because it is locally conservative 
to machine precision.) The error was calculated as div[i;] over each element. Notice that the error 
is concentrated along the inflow and outflow boundaries. This is due to the well-known implications 
of imposing the pressure boundary conditions in weak fashion. 



Remark 3. Due to the well-known instabilities associated with a high ratio of advection to diffusion 
in the transport model |321 133] . some diffusion was necessary to include for numerical reasons, even 
though the primary quantities of interest are related to pure advection. Also a constant scalar 
diffusivity was chosen to avoid instabilities associated with violations of the discrete maximum- 
minimum principle [34^ \35\ l36] . 
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Table 2. Multilayer problem: parameters 



a T* a TTI PT" 


\/ a 1 n p 

V CXi U.C 




QQfi 1 Vcr/m^ 


ft/1 


J. 111 




5x10^^^ ir? 


"•A 


n 5x10^^^ 




Sxin^^^ 




o 1 n— 1 2 

8x10 ^■^ m^ 




1x10^^ kg/m-s 




(0,0,0) 




1x10^ Pa 


Nodes 


441 


Elements 


800 



4.2. Cylinder inclusion problem. The cylinder inclusion problem can be conceptualized as a 
contaminant initially contained in a highly permeable porous media, encapsulated in a lower per- 
meability media, that is dispersed by a cross-flow in the domain generated by a pressure drop. This 
problem is an adaptation from the problem presented in j37] . The primary quantity of interest 
in this problem is the degree to which the contaminant infiltrates the domain and the time scale 
of this process. The geometry of this problem is shown in Figure |8] along with the computational 
mesh. The permeability of the infinite domain, and other parameters is given Table [3j 

Figure [9] shows the velocity magnitude for both elements. For the RTO element, a smooth 
transition occurs from the highly permeable region to the region of lower permeability. On the 
other hand, the VMS element shows oscillations along the boundary of the two regions. Again, 
this is due to the continuity of the solution in the direction transverse to the interface for the VMS 



element. Figures 10 and 11 show the pressure and velocity along different portions of the domain. 
The overshoot of the velocity magnitude for the VMS element is also evident in these figures. 



Figure 12 shows the concentration of the species in time as it deviates from the initial conditions. 
Both the RTO and VMS elements again show tremendous similarity in the contours of concentration 
with respect to time. The local mass balance error for the VMS formulation is shown in Figure 



13 Unlike in the multilayer problem, for the cylinder inclusion, the error is concentrated along 



the boundary of the inclusion, in particular where the interface is oriented tangentially to the fiow 
direction. To more fully understand if these features are related to the magnitude of the jump 
in permeability, several cases were run in which the permeability of the inclusion was varied from 



1x10"^^ to 1x10-3 m^. Fi gure 14 shows the error in the concentration for both elements as the 
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jump increases. As expected, the RTO element performs well regardless of the jump in permeability, 
but even the VMS shows only a marginal error (less than 2%) even for the worst case scenario. 

Table 3. Cylinder inclusion problem: parameters 



Parameter 


Value 


P 


996.1 kg/m^ 


h 


1x10"^ m^ 




varies 1x10^^^ to lxlO~^ m^ 




1.13x10"^ kg/m-s 


9 


(0,0,0) 


Po 


2.1721x10^ Pa 


Nodes 


489 


Elements 


912 



4.3. Leaky well problem. The leaky well problem represents the transport of a species as it is 
injected into an aquifer that contains an abandoned well that allows some of the species to escape. 
The primary quantity of interest is the rate at which the species escapes from the aquifer through 
the leak. This problem has been studied by a number of authors in evaluating subsurface modeling 
codes (see [MIEHIIH]) and is therefore important to study from the perspective of element choice, 
since the element choice could have a significant effect on the results. The domain is shown in 



Figure 15 along with the computational mesh. 



The steady state pressure and velocity contours are shown in Figures 16 and [17) Given the 
velocity field computed by either element, the transport problem is solved by prescribing a Dirichlet 
boundary condition of c = 1.0 at the inflow. Thus, the injection rate is J cvy dAj and the leak 
rate is j cVy dAo, where Aj is the area of the injection port and Aq is the area of the leaky well, 
both with diameter 0.3 m. For the leak rate, c is measured at the opening of the leaky well into 



Aquifer A. Figures 18 and 19 show similar results for the overall contours of concentration for 



results obtained using the Darcy flow model given in equation (la). 

Figure 20 shows the leak rate for various values of /3 from 0.0 to 1x10"^ using the Barus model 



given in equation ( lb ) . Figure 20 shows that the dependence of the viscosity on the pressure has a 



strong influence on the predicted leak rate. This trend is also apparent in Figure 21 , which shows 
contours of the concentration at time t = 0.12 days for various values of (3. Notice the greater 
degree of penetration for the low values of /3 vs. high. The pressure dependent viscosity leads to 



both a lower leak rate and a longer time to steady-state conditions. Figure 20 reveals that the 
leak rate is substantially under-predicted by the VMS element. Although both elements predict 
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that the leak rate will reach a steady-state around t = 2 days for the case of /? = 0.0, the VMS 
element shows a considerable decrease in the steady-state leak rate. Whereas for the previous two 
problems, the RTO and VMS element performed remarkably similar, in this case the VMS element 
shows considerable variation. The results suggest that the VMS element is not an appropriate 
choice for this type of problem since it not only shows a significant accuracy error, but also an 
under-prediction of the quantity of interest which, from a design standpoint, is less conservative. 

Table 4. Leaky well problem: parameters 



Parameter Value 



p 479 kg/m^ 

ki 1x10^^2 

k2 IxlO-i-^m^ 

IJ,0 3.95x10"^ kg/m-s 

g (0,0,-9.8) m/s^ 

Injection pressure 2.03x10^ Pa 

Open boundary pressure 3.075x10^ + 1.025x10^ (depth) Pa 

Nodes 23126 

Elements 43098 



4.4. Computational efficiency. In this section, we evaluate both elements from the perspective 

of computation time and memory usage. The most significant drawback pointed out for the RTO 

element is that it is not scalable as the problem size grows. This is due to the RTO element 

formulation having variables that exist on the element edges, which in turn leads to more equations 

in the linear system. Simply put, a node unknown can be shared by many elements, where as 

an edge or face unknown can be shared by only two elements. The fewer elements that can 

share an unknown, the larger the system. Although the problems shown in this work are only two 

dimensional, effects of this drawback are present. Table[5]shows the number of unknowns, assembly 

time, solve time, and memory usage of both elements for all of the examples above. Note that for 

each problem the number of unknowns is substantially higher for the RTO element, roughly 30%. 

If this were the only consideration, the VMS element would be the clear winner, but taking into 

account the computation time as well, the VMS element takes substantially longer to assemble and 

solve the linear system. This is due to the additional stabilization terms that must be computed 

at each step. As expected, the time savings decrease as the problem size increases since at some 

point the problem size offsets the faster computation for the RTO element. 
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The same is true for the memory usage. For the two smaller problems, the memory usage is 
identical, but for the leaky well problem, the VMS element requires less memory. The results 
suggest that for small to moderately sized problems (less than 100,000 degrees of freedom) the 
RTO element is more efficient, but for large problems the VMS element is more efficient. 

Table 5. Timing and memory usage data 

Multilayer Cylinder Inclusion Leaky well 



Quantity RTO VMS RTO VMS RTO VMS 



Total unknowns 


2481 


1764 


2801 


1956 


132447 


92504 


Assembly time/time step (s) 


0.0323 


0.05162 


0.0345 


0.05746 


2.943 


4.246 


System solve/time step (s) 


0.0327 


0.04056 


0.0156 


0.04959 


2.564 


2.320 


Required memory storage (MB) 


233 


233 


227 


227 


421 


396 


Problem size savings 




28.9% 




30.1% 




30.2% 


Time savings 


29.5% 




53.2% 




16.1% 




Memory savings 












5.9% 



5. CONCLUSIONS 

The results of this study reveal several interesting features of the lowest Raviart-Thomas and the 
variational multiscale formulation when used to model flow through heterogeneous porous media. 
Although the equal-order interpolation for the velocity and pressure under the variational multiscale 
formulation does not exhibit local (or element) mass balance. In addition, the formulation also 
shows oscillations in the velocity field when the primary flow direction is tangential to jumps in 
permeability, for some transport problems, the resulting concentration field is only marginally 
affected. For the multilayer problem and the cylinder inclusion problem, the mass conservation 
error produced by using the variational multiscale formulation is less than 2%. On the other hand, 
for transport problems similar to the leaky well problem, the equal-order interpolation under the 
variational multiscale formulation not only shows a substantial error in the concentration, but 
also under-predicts the leak rate, which is undesirable from a design perspective. For the leaky 
well problem, the leak rate was under-predicted by roTigiily 10% by the variational multiscale 
formulation. The results of this study also reveal that for small to moderately sized problems, 
the increased size of the linear system of equations for the RTO formulation is offset by a reduced 
computation time. For large problems, this advantage diminishes. In general, these results give 
some perspective on the computational cost associated with using a locally mass conserving element, 
and the accuracy cost of using a stabilized element. 
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V ■ n 

Figure 1. The left figure shows the degrees-of- freedom in the lowest-order Raviart- 
Thomas (RTO) formulation, and the right figure shows the degrees-of-freedom in the 
(equal-order) VMS formulation. 




Figure 2. Multilayer problem: The left figure shows the geometry and boundary 
conditions, and the right figure shows the computational mesh employed in the 
numerical simulation. 
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Figure 3. Multilayer problem: Variation of x- velocity along the vertical centerline 
of the domain at x = 0.5. 
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Figure 4. Multilayer problem: Variation of the concentration along the vertical 
line of the domain at x = 1.0 for various instances of time t = 0.5, 1.0 and 1.5 s. 
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Figure 5. Multilayer problem: Contours of the concentration at time t = 0.2 s 
using the RTO formulation (left) and the variational multiscale formulation (right). 
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Figure 6. Multilayer problem: This figure shows the variation of the total concen- 
tration with time. The total concentration is the integral of the concentration field 
over the entire domain. 
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Figure 7. Multilayer problem: This figure shows the contours of the error in the 
local mass balance under the equal-order variational multiscale formulation. 
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Figure 8. Cylinder inclusion: The top figure shows the problem geometry and 
boundary conditions, and the bottom figure shows the computational mesh used in 
the numerical simulation. 
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Figure 9. Cylinder inclusion: This figure shows the contours of the magnitude of 
the velocity using the RTO formulation (top) and the variational multiscale formu- 
lation (bottom) for the case k2 = 1 x 10^^. 
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Figure 10. Cylinder inclusion: This figure shows the variation of the pressure along 
x-direction at y = for the case /c2 = 1 x 10~^. 




Figure 11. Cylinder inclusion: This figure shows the variation of the x- velocity 
along the x-direction at y = (left), and along the y-direction at x = (right) for 
the case k2 = 1 x 10^^. 
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Figure 12. Cylinder inclusion: This figure shows the contours of the concentration 
at time t = 1.5 s using the RTO formulation (top) and the variational multiscale 
formulation (bottom) for the case k2 = I >^ 10^^. 
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Figure 13. Cylinder inclusion: This figure shows the contours of the error in the 
mass balance using the variational multiscale formulation for the case k2 = 1 x 10^^. 
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Figure 14. Cylinder inclusion: The figure shows the variation of the error in the 
total concentration (which is measured as the difference between the integral of c 
over the domain at time t = 2.5 s minus the integral of c over the domain at t = 0.0 s) 
with respect to the log of the ratio of permeabilities between the inclusion and the 
rest of the domain. 
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Figure 15. Leaky well: A pictorial description of the problem (top) and the com- 
putational mesh used in the numerical simulation (bottom). (Note that a portion 
of the mesh for the abandoned well is omitted to facilitate plotting.) 
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Figure 16. Leaky well: velocity magnitude contours using the RTO formulation 
(top) and the variational multiscale formulation (bottom). 
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Figure 17. Leaky well: pressure contours using the RTO formulation (top) and the 
variational multiscale formulation (bottom). 
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Figure 18. Leaky well: contours of concentration at various times for the RTO 
formulation. (Note that a portion of the abandoned well was omitted to facilitate 
plotting.) 
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Figure 19. Leaky well: contours of concentration at various instances of time for 
the variational multiscale formulation. (Note that a portion of the abandoned well 
was omitted to facilitate plotting.) 
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Figure 20. Leaky well: leak rate vs. time measured as the concentration times 
the mass flux across the opening of the abandoned well into Aquifer B for various 



values of /3 using the Barus model for viscosity given in equation (lb). 
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Figure 21. Leaky well: contours of the concentration at the injection well at t = 
0.12 days for various values of /3 using the Barus model for viscosity given in equation 



(lb) (top) RTO (bottom) VMS. 
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